# 1.DESeq差异分析
setwd("C:\\Users\\lenovo\\Desktop\\R语言")    #设置文件路径
data=read.csv("mRNA_exprSet.csv")        #导入数据
row.names(data)=data[,1]               #(以下两行)将基因名作为列名
data=data[-1]
data1 <- data[apply(data, 1, sum) > 0 , ]     #过滤全零数据
if (!require("BiocManager", quietly = TRUE))     #安装DESeq2包
+   install.packages("BiocManager")
library(DESeq2)
sample=colnames(data1)                #提取列名
design1=substring(sample,14,15)       #提取分组信息"01""11",设置design
coldata1=data.frame(sample,design1)    #设置coldata
dds <- DESeqDataSetFromMatrix(         #设置参数
  countData = data1,
  colData =coldata1,
  design = ~ design1)
dds$design1 <- relevel(dds$design1,ref="11")  #设置对照(正常)组
dds <- DESeq(dds)               #进行差异分析
res <- results(dds)             #输出差异分析结果
res <- cbind(sample1="tumour",sample2="normal",as.data.frame(res))   
res1=res[is.na(res$padj)==FALSE,]    #去除P.adj值为NA的基因
res2=res[is.na(res$pvalue)==FALSE,]     #去除Pvalue值为NA的基因 
res005 <- res1[res1$padj<0.05,]#提取padj<0.05的基因，以0.05为阈值，这些基因的P值可信度高
write.csv(res005,"res005.csv")

bgdata=res2                #此处以res2作为原始数据，绘制火山图
library("ggrepel")
library("ggplot2")
pthreshold=0.05             #(以下两行)设置P和FC的阈值
fcthreshold=1.2
# 设置点的分类
bgdata$change <- as.factor(ifelse(bgdata$pvalue<pthreshold & abs(bgdata$log2FoldChange)>log2(fcthreshold),
                                  ifelse(bgdata$log2FoldChange>log2(fcthreshold),"Up","Down"),"Non"))

# 样本标签
bgdata$label <- ifelse(bgdata$pvalue<pthreshold & abs(bgdata$log2FoldChange)>log2(fcthreshold),row.names(bgdata),"")

# 绘制火山图
p.vol <- ggplot(data = bgdata,
                aes(x =log2FoldChange ,y = -log10(pvalue),colour = change,fill = change))+
  scale_color_manual(values = c('green','grey','red'))+
  geom_point(alpha = 0.4,size = 3.5)+
  # 标签
  geom_text_repel(aes(x = log2FoldChange,y = -log10(pvalue),label = label),size = 3,
                  box.padding = unit(0.6,"lines"),point.padding = unit(0.7,"lines"),
                  segment.color = "black",show.legend = FALSE)+
  # 辅助线
  geom_vline(xintercept = c(-log2(fcthreshold),log2(fcthreshold)),lty = 4,col = "black",lwd = 0.8)+
  geom_hline(yintercept = -log10(pthreshold),lty = 4,col = "black",lwd = 0.8)+
  theme_bw()+
  labs(x = "log2(FoldChange)",y = "-log10(p)",title = "Volcano Plot of  Different Expression Proteins")+
  # 坐标轴标题、标签和图例相关设置
  theme(axis.text = element_text(size = 11),axis.title = element_text(size = 13), # 坐标轴标签和标题
        plot.title = element_text(hjust = 0.5,size = 15,face = "bold"), # 标题
        legend.text = element_text(size = 11),legend.title = element_text(size = 13), # 图例标签和标题
        plot.margin = unit(c(0.5,0.5,0.5,0.5),"cm")) # 图边距
ggsave(p.vol,filename = "Gene Diference analysis1.pdf")
